Using the Benford’s law to assess the quality of COVID-19 reported data in Vietnam
Preambule
Analysis following previously published work on Brazilian data, using the benford.analysis R package.
Conclusions
To be completed.
Packages and global variables
clean_folder <- "~/Dropbox/aaa/studies/oucomo/clean data/"library(readxl)
library(magrittr)
library(dplyr)
library(tidyr)
library(purrr)
library(benford.analysis)
library(tibble)
library(lubridate)
library(sf)
library(stringr)
library(stringi)Utilitary functions
A function to extract the statistics of a Benford analysis we want to look at:
extract_statistics <- function(x) {
with(x, c(with(info, c(n, n.second.order)),
stats$chisq$statistic,
MAD,
MAD.conformity,
distortion.factor,
stats$mantissa.arc.test$statistic))
}A function to make a table of statistics (in column) for each of the provinces (in row):
make_statistics_table <- function(x) {
x %>%
map_df(extract_statistics) %>%
t() %>%
as.data.frame() %>%
rownames_to_column() %>%
setNames(c("province", "n", "n2", "Chisq", "MAD", "MAD conformity", "DF", "Mantissa")) %>%
as_tibble() %>%
mutate_at(vars(starts_with("n")), as.integer) %>%
mutate_at(vars(MAD, DF, Chisq, Mantissa), as.numeric)
}A function that plots the diagnostic plots for each of the provinces:
plot_digits_distribution <- function(x, y) {
xmarks <- barplot(x[["bfd"]]$data.dist.freq, col = 4,
xlab = "digits", ylab = "frequency",
ylim = c(0, 1.1 * max(c(x[["bfd"]]$data.dist.freq,
x[["bfd"]]$benford.dist.freq))))
axis(1, at = xmarks, labels = x[["bfd"]]$digits)
lines(xmarks, x[["bfd"]]$benford.dist.freq, lwd = 2, col = 2)
title(y, line = -1)
}Data
Vietnam 2019 census data
The 2019 census data (Vinh Long is missing):
census <- paste0(clean_folder, "census2019.rds") %>%
readRDS() %>%
group_by(province) %>%
summarise(popsize = sum(n)) %>%
mutate_at("province", stri_trans_general, "Latin-ASCII") %>%
mutate_at("province", str_remove, "Thanh pho |Tinh ") %>%
bind_rows(tibble(province = "Vinh Long", popsize = 1141677))Geographic data of Vietnam
Downloading the file if not in the folder:
if(!file.exists("gadm36_VNM_1_sf.rds"))
download.file("https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/gadm36_VNM_1_sf.rds",
"gadm36_VNM_1_sf.rds")Loading the data:
gadm <- "gadm36_VNM_1_sf.rds" %>%
readRDS() %>%
st_set_crs(4326) %>%
transmute(province = VARNAME_1, geometry) %>%
mutate(centroid = map(geometry, st_centroid),
coordinates = map(centroid, st_coordinates),
latitude = map(coordinates, extract2, 2)) %>%
left_join(census, "province")Rogier’s COVID-19 data
Reading the data:
rogier <- paste0(clean_folder, "Fourth cluster First wave.xlsx") %>%
read_excel() %>%
rename(date = `...1`)Getting rid off whatever is below the Total row:
rogier %<>% head(which(rogier$date == "Total") - 1)Reformatting the data:
rogier %<>%
select(date, `An Giang`:last_col()) %>%
filter(! is.na(date)) %>%
mutate_all(replace_na, 0) %>%
mutate_all(as.integer) %>%
mutate_at("date", as.Date, origin = "1899-12-30")which looks like:
rogier## # A tibble: 339 x 64
## date `An Giang` `Ba Ria - Vung Tau` `Bac Giang` `Bac Kan` `Bac Lieu`
## <date> <int> <int> <int> <int> <int>
## 1 2020-01-23 0 0 0 0 0
## 2 2020-02-01 0 0 0 0 0
## 3 2020-02-03 0 0 0 0 0
## 4 2020-02-04 0 0 0 0 0
## 5 2020-02-06 0 0 0 0 0
## 6 2020-02-07 0 0 0 0 0
## 7 2020-02-09 0 0 0 0 0
## 8 2020-02-11 0 0 0 0 0
## 9 2020-02-13 0 0 0 0 0
## 10 2020-03-06 0 0 0 0 0
## # … with 329 more rows, and 58 more variables: Bac Ninh <int>, Ben Tre <int>,
## # Binh Dinh <int>, Binh Duong <int>, Binh Phuoc <int>, Binh Thuan <int>,
## # Ca Mau <int>, Can Tho <int>, Cao Bang <int>, Da Nang <int>, Dak Lak <int>,
## # Dak Nong <int>, Dien Bien <int>, Dong Nai <int>, Dong Thap <int>,
## # Gia Lai <int>, Ha Giang <int>, Ha Nam <int>, Ha Noi <int>, Ha Tinh <int>,
## # Hai Duong <int>, Hai Phong <int>, Hau Giang <int>, Ho Chi Minh <int>,
## # Hoa Binh <int>, Hue <int>, Hung Yen <int>, Khanh Hoa <int>,
## # Kien Giang <int>, Kon Tum <int>, Lai Chau <int>, Lam Dong <int>,
## # Lang Son <int>, Lao Cai <int>, Long An <int>, Nam Dinh <int>,
## # Nghe An <int>, Ninh Binh <int>, Ninh Thuan <int>, Phu Tho <int>,
## # Phu Yen <int>, Quang Binh <int>, Quang Nam <int>, Quang Ngai <int>,
## # Quang Ninh <int>, Quang Tri <int>, Soc Trang <int>, Son La <int>,
## # Tay Ninh <int>, Thai Binh <int>, Thai Nguyen <int>, Thanh Hoa <int>,
## # Tien Giang <int>, Tra Vinh <int>, Tuyen Quang <int>, Vinh Long <int>,
## # Vinh Phuc <int>, Yen Bai <int>For compatibility of other data sets, we want to rename Hue into Thua Thien Hue:
names(rogier) <- str_replace(names(rogier), "Hue", "Thua Thien Hue")NCSC’s COVID-19 data
To be completed.
Analyzing Rogier’s data
All the data
Removing the dates and filtering out the zeros:
rogier1 <- rogier %>%
select(-date) %>%
map(~ .x[.x > 0])Performing the Benford analysis on each of the 63 provinces:
benford_analyses_rogier1 <- map(rogier1, benford, number.of.digits = 1)Extracting all these statistics for each of the 63 provinces. Interpretation is as follows:
n: number of data points on which the statistics are calculated.n2: number of data points on which the second order test is performed.Chisq: classical chi-square statistic. Note though that it is highly sensitive to the sample size and tends to reject the null even for small departures from the expected distribution.MADstatistic is thus to be preferred.MAD: Mean Absolute Deviation. The MAD test is more robust since it ignores the number of records. The higher the MAD, the larger the average difference between the observed and theoretical distributions. MAD values above 0.015 suggest nonconformity.MAD conformity: interpretation of theMADstatistic value.DF: Distortion Factor. The DF statistic examines the digit patterns to indicate whether the data appear to be underestimated or overestimated and the deformity’s magnitude.Mantissa: expected to be 0.5. Values less than 0.5 suggest that figures tend to be manipulated downward whereas it’s the opposite for values higher than 0.5.
table_rogier1 <- make_statistics_table(benford_analyses_rogier1)
print(table_rogier1, n = Inf)## # A tibble: 63 x 8
## province n n2 Chisq MAD `MAD conformity` DF Mantissa
## <chr> <int> <int> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 An Giang 143 114 4.38 0.0145 Marginally acceptabl… -13.7 0.00177
## 2 Ba Ria - V… 129 82 16.1 0.0342 Nonconformity -5.98 0.00984
## 3 Bac Giang 123 49 10.5 0.0243 Nonconformity -15.0 0.00331
## 4 Bac Kan 13 2 12.6 0.0996 Nonconformity NaN 0.274
## 5 Bac Lieu 124 68 19.4 0.0373 Nonconformity -16.5 0.0390
## 6 Bac Ninh 139 54 1.29 0.00899 Acceptable conformity -15.3 0.00597
## 7 Ben Tre 125 73 9.08 0.0238 Nonconformity -5.88 0.00707
## 8 Binh Dinh 132 58 17.7 0.0341 Nonconformity -37.5 0.0676
## 9 Binh Duong 147 142 60.0 0.0550 Nonconformity 21.6 0.135
## 10 Binh Phuoc 115 66 10.8 0.0297 Nonconformity -29.9 0.0176
## 11 Binh Thuan 121 86 44.0 0.0543 Nonconformity 7.89 0.139
## 12 Ca Mau 112 68 12.3 0.0336 Nonconformity -9.55 0.0392
## 13 Can Tho 133 90 22.0 0.0364 Nonconformity 9.13 0.00876
## 14 Cao Bang 23 15 8.82 0.0629 Nonconformity -65.1 0.0746
## 15 Da Nang 120 64 7.92 0.0205 Nonconformity -0.793 0.0221
## 16 Dak Lak 106 80 12.2 0.0329 Nonconformity -11.8 0.0500
## 17 Dak Nong 115 45 14.1 0.0312 Nonconformity -4.75 0.0584
## 18 Dien Bien 45 19 5.21 0.0322 Nonconformity -44.4 0.0227
## 19 Dong Nai 151 134 101. 0.0808 Nonconformity 50.6 0.277
## 20 Dong Thap 142 102 35.5 0.0423 Nonconformity 8.72 0.0367
## 21 Gia Lai 109 57 7.18 0.0203 Nonconformity -1.53 0.00845
## 22 Ha Giang 51 45 28.2 0.0797 Nonconformity -6.27 0.202
## 23 Ha Nam 90 31 30.5 0.0594 Nonconformity -45.2 0.104
## 24 Ha Noi 168 78 3.56 0.0127 Marginally acceptabl… -2.41 0.0172
## 25 Ha Tinh 99 29 17.1 0.0409 Nonconformity -36.2 0.00290
## 26 Hai Duong 112 36 15.2 0.0322 Nonconformity -25.8 0.00757
## 27 Hai Phong 45 17 11.2 0.0489 Nonconformity -48.8 0.0412
## 28 Hau Giang 104 62 5.61 0.0190 Nonconformity -19.4 0.00799
## 29 Ho Chi Minh 192 174 23.5 0.0308 Nonconformity -5.19 0.00437
## 30 Hoa Binh 32 15 12.8 0.0520 Nonconformity 5.36 0.000277
## 31 Thua Thien… 103 51 11.3 0.0299 Nonconformity -15.2 0.0444
## 32 Hung Yen 65 23 19.9 0.0547 Nonconformity -42.9 0.112
## 33 Khanh Hoa 137 90 16.8 0.0286 Nonconformity -9.27 0.0351
## 34 Kien Giang 137 103 30.1 0.0392 Nonconformity 6.00 0.0106
## 35 Kon Tum 62 18 6.97 0.0298 Nonconformity -60.5 0.00649
## 36 Lai Chau 16 4 11.4 0.0755 Nonconformity NaN 0.146
## 37 Lam Dong 92 42 8.39 0.0244 Nonconformity -14.7 0.00723
## 38 Lang Son 77 16 3.86 0.0223 Nonconformity -44.8 0.0370
## 39 Lao Cai 49 7 8.75 0.0319 Nonconformity -54.0 0.00321
## 40 Long An 144 113 17.4 0.0282 Nonconformity 12.8 0.0564
## 41 Nam Dinh 69 40 10.7 0.0401 Nonconformity -3.09 0.0493
## 42 Nghe An 124 63 6.02 0.0180 Nonconformity -4.56 0.00912
## 43 Ninh Binh 46 12 7.86 0.0353 Nonconformity -61.7 0.0185
## 44 Ninh Thuan 117 48 20.3 0.0350 Nonconformity -10.4 0.0776
## 45 Phu Tho 65 37 16.6 0.0435 Nonconformity 6.64 0.0671
## 46 Phu Yen 153 47 10.9 0.0260 Nonconformity -37.7 0.0181
## 47 Quang Binh 104 49 9.00 0.0242 Nonconformity -28.6 0.00553
## 48 Quang Nam 108 51 8.04 0.0237 Nonconformity -6.50 0.0171
## 49 Quang Ngai 119 47 15.9 0.0285 Nonconformity -26.9 0.00436
## 50 Quang Ninh 54 27 13.2 0.0432 Nonconformity -32.6 0.0500
## 51 Quang Tri 88 27 3.57 0.0164 Nonconformity -48.3 0.0144
## 52 Soc Trang 82 67 6.18 0.0261 Nonconformity -11.6 0.000278
## 53 Son La 66 18 15.0 0.0454 Nonconformity -60.9 0.0710
## 54 Tay Ninh 131 110 20.4 0.0307 Nonconformity 6.56 0.0342
## 55 Thai Binh 62 26 5.63 0.0248 Nonconformity 20.6 0.0784
## 56 Thai Nguyen 38 16 5.65 0.0306 Nonconformity 24.2 0.0156
## 57 Thanh Hoa 117 48 9.67 0.0291 Nonconformity -3.92 0.0191
## 58 Tien Giang 139 118 8.39 0.0229 Nonconformity -5.77 0.00817
## 59 Tra Vinh 119 76 9.03 0.0239 Nonconformity -9.04 0.00150
## 60 Tuyen Quang 40 11 54.4 0.0942 Nonconformity -26.6 0.178
## 61 Vinh Long 137 74 24.1 0.0307 Nonconformity -8.72 0.0264
## 62 Vinh Phuc 60 31 8.82 0.0378 Nonconformity -22.6 0.0265
## 63 Yen Bai 28 8 18.5 0.0844 Nonconformity -49.9 0.206Plotting the diagnostic plots for each of the 63 provinces:
par(mfrow = c(16, 4), plt = c(.2, .97, .17, .9))
walk2(benford_analyses_rogier1,
names(benford_analyses_rogier1),
plot_digits_distribution)
Data from January 2021
Same analysis as above but with data from January 2021 only:
rogier2 <- rogier %>%
filter(date > ymd(20210101)) %>%
select(-date) %>%
map(~ .x[.x > 0])
benford_analyses_rogier2 <- map(rogier2, benford, number.of.digits = 1)
table_rogier2 <- make_statistics_table(benford_analyses_rogier2)
print(table_rogier2, n = Inf)## # A tibble: 63 x 8
## province n n2 Chisq MAD `MAD conformity` DF Mantissa
## <chr> <int> <int> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 An Giang 143 114 4.38 0.0145 Marginally acceptabl… -13.7 0.00177
## 2 Ba Ria - V… 129 82 16.1 0.0342 Nonconformity -5.98 0.00984
## 3 Bac Giang 123 49 10.5 0.0243 Nonconformity -15.0 0.00331
## 4 Bac Kan 13 2 12.6 0.0996 Nonconformity NaN 0.274
## 5 Bac Lieu 124 68 19.4 0.0373 Nonconformity -16.5 0.0390
## 6 Bac Ninh 139 54 1.29 0.00899 Acceptable conformity -15.3 0.00597
## 7 Ben Tre 125 73 9.08 0.0238 Nonconformity -5.88 0.00707
## 8 Binh Dinh 132 58 17.7 0.0341 Nonconformity -37.5 0.0676
## 9 Binh Duong 147 142 60.0 0.0550 Nonconformity 21.6 0.135
## 10 Binh Phuoc 115 66 10.8 0.0297 Nonconformity -29.9 0.0176
## 11 Binh Thuan 121 86 44.0 0.0543 Nonconformity 7.89 0.139
## 12 Ca Mau 112 68 12.3 0.0336 Nonconformity -9.55 0.0392
## 13 Can Tho 133 90 22.0 0.0364 Nonconformity 9.13 0.00876
## 14 Cao Bang 23 15 8.82 0.0629 Nonconformity -65.1 0.0746
## 15 Da Nang 119 64 7.73 0.0203 Nonconformity -0.793 0.0203
## 16 Dak Lak 106 80 12.2 0.0329 Nonconformity -11.8 0.0500
## 17 Dak Nong 115 45 14.1 0.0312 Nonconformity -4.75 0.0584
## 18 Dien Bien 45 19 5.21 0.0322 Nonconformity -44.4 0.0227
## 19 Dong Nai 151 134 101. 0.0808 Nonconformity 50.6 0.277
## 20 Dong Thap 142 102 35.5 0.0423 Nonconformity 8.72 0.0367
## 21 Gia Lai 109 57 7.18 0.0203 Nonconformity -1.53 0.00845
## 22 Ha Giang 51 45 28.2 0.0797 Nonconformity -6.27 0.202
## 23 Ha Nam 90 31 30.5 0.0594 Nonconformity -45.2 0.104
## 24 Ha Noi 165 78 3.32 0.0117 Acceptable conformity -2.41 0.0165
## 25 Ha Tinh 99 29 17.1 0.0409 Nonconformity -36.2 0.00290
## 26 Hai Duong 112 36 15.2 0.0322 Nonconformity -25.8 0.00757
## 27 Hai Phong 45 17 11.2 0.0489 Nonconformity -48.8 0.0412
## 28 Hau Giang 104 62 5.61 0.0190 Nonconformity -19.4 0.00799
## 29 Ho Chi Minh 188 174 23.6 0.0316 Nonconformity -5.19 0.00345
## 30 Hoa Binh 32 15 12.8 0.0520 Nonconformity 5.36 0.000277
## 31 Thua Thien… 103 51 11.3 0.0299 Nonconformity -15.2 0.0444
## 32 Hung Yen 65 23 19.9 0.0547 Nonconformity -42.9 0.112
## 33 Khanh Hoa 137 90 16.8 0.0286 Nonconformity -9.27 0.0351
## 34 Kien Giang 137 103 30.1 0.0392 Nonconformity 6.00 0.0106
## 35 Kon Tum 62 18 6.97 0.0298 Nonconformity -60.5 0.00649
## 36 Lai Chau 16 4 11.4 0.0755 Nonconformity NaN 0.146
## 37 Lam Dong 92 42 8.39 0.0244 Nonconformity -14.7 0.00723
## 38 Lang Son 77 16 3.86 0.0223 Nonconformity -44.8 0.0370
## 39 Lao Cai 49 7 8.75 0.0319 Nonconformity -54.0 0.00321
## 40 Long An 144 113 17.4 0.0282 Nonconformity 12.8 0.0564
## 41 Nam Dinh 69 40 10.7 0.0401 Nonconformity -3.09 0.0493
## 42 Nghe An 124 63 6.02 0.0180 Nonconformity -4.56 0.00912
## 43 Ninh Binh 46 12 7.86 0.0353 Nonconformity -61.7 0.0185
## 44 Ninh Thuan 117 48 20.3 0.0350 Nonconformity -10.4 0.0776
## 45 Phu Tho 65 37 16.6 0.0435 Nonconformity 6.64 0.0671
## 46 Phu Yen 153 47 10.9 0.0260 Nonconformity -37.7 0.0181
## 47 Quang Binh 104 49 9.00 0.0242 Nonconformity -28.6 0.00553
## 48 Quang Nam 108 51 8.04 0.0237 Nonconformity -6.50 0.0171
## 49 Quang Ngai 119 47 15.9 0.0285 Nonconformity -26.9 0.00436
## 50 Quang Ninh 54 27 13.2 0.0432 Nonconformity -32.6 0.0500
## 51 Quang Tri 88 27 3.57 0.0164 Nonconformity -48.3 0.0144
## 52 Soc Trang 82 67 6.18 0.0261 Nonconformity -11.6 0.000278
## 53 Son La 66 18 15.0 0.0454 Nonconformity -60.9 0.0710
## 54 Tay Ninh 131 110 20.4 0.0307 Nonconformity 6.56 0.0342
## 55 Thai Binh 62 26 5.63 0.0248 Nonconformity 20.6 0.0784
## 56 Thai Nguyen 38 16 5.65 0.0306 Nonconformity 24.2 0.0156
## 57 Thanh Hoa 117 48 9.67 0.0291 Nonconformity -3.92 0.0191
## 58 Tien Giang 139 118 8.39 0.0229 Nonconformity -5.77 0.00817
## 59 Tra Vinh 119 76 9.03 0.0239 Nonconformity -9.04 0.00150
## 60 Tuyen Quang 40 11 54.4 0.0942 Nonconformity -26.6 0.178
## 61 Vinh Long 137 74 24.1 0.0307 Nonconformity -8.72 0.0264
## 62 Vinh Phuc 60 31 8.82 0.0378 Nonconformity -22.6 0.0265
## 63 Yen Bai 28 8 18.5 0.0844 Nonconformity -49.9 0.206par(mfrow = c(16, 4), plt = c(.2, .97, .17, .9))
walk2(benford_analyses_rogier2,
names(benford_analyses_rogier2),
plot_digits_distribution)
Let’s compare with the table we get with all the data:
vars <- setdiff(names(table_rogier1), c("province", "MAD conformity"))
par(mfrow = c(2, 3), plt = c(.2, .97, .17, .9), cex = 1)
walk(vars, ~ {
plot(table_rogier1[[.x]], table_rogier2[[.x]], col = 4,
xlab = "all the data", ylab = "from Jan 2021")
abline(0, 1)
title(.x, line = -1)})
Specific comparisons
Group 1: Hai Duong and Quang Ninh between 1 January 2021 and 1 April 2021
group1 <- rogier %>%
transmute(date, n = `Hai Duong` + `Quang Ninh`) %>%
filter(ymd(20210101) < date, date < ymd(20210401)) %>%
select(-date) %>%
map(~ .x[.x > 0])Group 2: Bac Giang and Bac Ninh between 1 April 2021 and 1 August 2021
group2 <- rogier %>%
transmute(date, n = `Bac Giang` + `Bac Ninh`) %>%
filter(ymd(20210401) < date, date < ymd(20210801)) %>%
select(-date) %>%
map(~ .x[.x > 0])Group 3: HCMC and Mekong between 1 May 2021 and 15 October 2021
The fourth wave:
wave4 <- rogier %>%
filter(ymd(20210501) < date, date < ymd(20211015)) %>%
select(-date)Merging with GADM data:
wave4gadm <- wave4 %>%
colSums() %>%
list() %>%
data.frame() %>%
setNames("n") %>%
rownames_to_column("province") %>%
left_join(gadm, ., "province") %>%
mutate(incidence = 100000 * n / popsize)The south wave in the south:
wave4mekong <- wave4gadm %>%
filter(latitude < 12, incidence > 50)Which looks like this (in red):
par(plt = c(0, 1, 0, 1))
gadm %>%
st_geometry() %>%
plot(col = "grey")
wave4mekong %>%
st_geometry() %>%
plot(add = TRUE, col = "red")
Group 3 is then:
group3 <- wave4 %>%
rowSums() %>%
list()Group 4: the whole country after 15 October 2021
group4 <- rogier %>%
filter(date > ymd(20211015)) %>%
select(-date) %>%
rowSums() %>%
list()Let’s group the 4 groups into a list:
groups <- c(group1, group2, group3, group4) %>%
setNames(paste0("group", 1:4))on which we can rerun the previous analyses:
benford_analyses_rogier_4groups <- map(groups, benford, number.of.digits = 1)
make_statistics_table(benford_analyses_rogier_4groups)## # A tibble: 4 x 8
## province n n2 Chisq MAD `MAD conformity` DF Mantissa
## <chr> <int> <int> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 group1 41 16 17.6 0.0465 Nonconformity -30.7 0.00933
## 2 group2 62 37 14.0 0.0481 Nonconformity -2.84 0.0852
## 3 group3 165 154 24.3 0.0282 Nonconformity 10.1 0.0329
## 4 group4 49 48 21.8 0.0563 Nonconformity 25.9 0.113par(mfrow = c(1, 4), plt = c(.2, .97, .17, .9))
walk2(benford_analyses_rogier_4groups,
names(benford_analyses_rogier_4groups),
plot_digits_distribution)
Analyzing NCSC’s data
To be done.